--- title: Peak slice maps (working on it) keywords: fastai sidebar: home_sidebar summary: "Slicing the cube " description: "Slicing the cube " nb_path: "notebooks/60_peakmaps.ipynb" ---
{% raw %}
{% endraw %} {% raw %}
import skimage.exposure as ske
{% endraw %} {% raw %}
from maxrf4u import DataStack, HotmaxAtlas, get_hotslices, get_peakmaps, multi_plot, plot_peakslices  

ds = DataStack('RP-T-1898-A-3689.datastack') 

x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum') 
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big) 
imvis_reg_highres = ds.read('imvis_reg_highres') 
imvis_reg = ds.read('imvis_reg') 
imvis_extent = ds.read('imvis_extent')

hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
{% endraw %}

How about Calcium?

{% raw %}
hma.hotmax_pixels[:,2]
array([  95,  151,  335,  359,  427,  465,  497,  547,  581,  685,  736,
        802,  900,  960,  984, 1017, 1150, 1261, 1358, 1539, 1573, 1976,
       2112, 2192])
{% endraw %} {% raw %}
def hotslice_indices(hma=hma): 
    '''Temporary hack to locate hotpixel slice index for  (indicated with a star in the multiplot).
    
    In the future do better ordering by putting hotslice first.'''
    
    keV_idxs = hma.hotmax_pixels[:,2]  
    
    hotslice_idxs = []
    
    for n, peak_idxs in enumerate(hma.peak_idxs_list): 
    
         hotslice_idxs.append(peak_idxs.index(keV_idxs[n])) 
    
    return hotslice_idxs 
{% endraw %} {% raw %}
hotslice_indices()[9]
4
{% endraw %} {% raw %}
def hotslice_index(n, hma=hma): 
    '''Temporary hack to locate hotpixel slice (indicated with a star).
    
    In the future do better ordering by putting hotslice first.'''
    
    keV_i = hma.hotmax_pixels[:,2][n]  
    
    hotslice_i = hma.peak_idxs_list[n].index(keV_i) 
    
    return hotslice_i
     
{% endraw %} {% raw %}
hotslice_index(3)
3
{% endraw %} {% raw %}
n = 5

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/9
[########################################] | 100% Completed |  6.5s
Computing slice 1/9
[########################################] | 100% Completed |  5.3s
Computing slice 2/9
[########################################] | 100% Completed | 12.7s
Computing slice 3/9
[########################################] | 100% Completed |  6.2s
Computing slice 4/9
[########################################] | 100% Completed |  5.5s
Computing slice 5/9
[########################################] | 100% Completed |  5.7s
Computing slice 6/9
[########################################] | 100% Completed |  6.4s
Computing slice 7/9
[########################################] | 100% Completed |  6.3s
Computing slice 8/9
[########################################] | 100% Completed |  5.7s
Computing slice 9/9
[########################################] | 100% Completed |  5.9s
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
titles[n]= '*' + titles[n]

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

fig, axs = multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel)
{% endraw %}

Let's check if peak #5[6] is an escape peak for Ca_Ka by calculating the energy shift....

{% raw %}
peak_idxs = hma.peak_idxs_list[5]
{% endraw %} {% raw %}
x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
1.7459311050874269
{% endraw %}

Indeed!!

{% raw %}
title = axs.flatten()[6].get_title() + '(escape)'
axs.flatten()[6].set_title(title)
Text(0.5, 1.0, '[6] 1.94 keV(escape)')
{% endraw %}

There is much more to be learned from this overview, but before going into this let's look at other elements.

How about Manganese?

Let's create a similar peak map overview for the hotmax spectrum corresponding to Fe_Ka. This is hotmax spectrum #10. Also of interest are the related hotmax spectra for manganese #9 and for Fe_Kb #11.

{% raw %}
n = 9 # Mn

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/10
[########################################] | 100% Completed |  5.4s
Computing slice 1/10
[########################################] | 100% Completed |  3.8s
Computing slice 2/10
[########################################] | 100% Completed |  6.4s
Computing slice 3/10
[########################################] | 100% Completed |  5.8s
Computing slice 4/10
[########################################] | 100% Completed |  6.1s
Computing slice 5/10
[########################################] | 100% Completed |  5.7s
Computing slice 6/10
[########################################] | 100% Completed |  6.2s
Computing slice 7/10
[########################################] | 100% Completed |  5.8s
Computing slice 8/10
[########################################] | 100% Completed |  6.6s
Computing slice 9/10
[########################################] | 100% Completed |  5.9s
Computing slice 10/10
[########################################] | 100% Completed |  6.1s
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

Mm, these histogram equalized peak maps aren't very informative. Most of the slices are just noise. I can not be sure how manganese and iron are correlated. Need to plot that in a different way.

{% raw %}
FeKa_map = peak_maps[2]
MnKa_map = peak_maps[4]
{% endraw %} {% raw %}
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])

ax.imshow(MnKa_map)
ax.set_title('Mn_Ka')

y, x, z = hot_pixel
ax.scatter([x], [y], c='r')
ax1.imshow(FeKa_map, vmax=2)
ax1.set_title('Fe_Ka')
ax1.scatter([x], [y], c='r');
{% endraw %}

My first question is, are there any other manganese speckles? And second, do the correlate with iron speckles, as I would expect?

This question needs to wait. Need speckle segmentation functions.

How about Iron?

{% raw %}
n = 10 # Fe_Ka

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/6
[########################################] | 100% Completed |  4.8s
Computing slice 1/6
[########################################] | 100% Completed |  6.0s
Computing slice 2/6
[########################################] | 100% Completed |  5.3s
Computing slice 3/6
[########################################] | 100% Completed |  6.2s
Computing slice 4/6
[########################################] | 100% Completed |  5.7s
Computing slice 5/6
[########################################] | 100% Completed |  6.5s
Computing slice 6/6
[########################################] | 100% Completed |  6.4s
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

How about Chlorine?

{% raw %}
n = 3

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/10
[########################################] | 100% Completed |  5.4s
Computing slice 1/10
[########################################] | 100% Completed | 11.6s
Computing slice 2/10
[########################################] | 100% Completed |  6.5s
Computing slice 3/10
[########################################] | 100% Completed |  3.4s
Computing slice 4/10
[########################################] | 100% Completed |  6.1s
Computing slice 5/10
[########################################] | 100% Completed |  6.3s
Computing slice 6/10
[########################################] | 100% Completed |  8.1s
Computing slice 7/10
[########################################] | 100% Completed | 11.8s
Computing slice 8/10
[########################################] | 100% Completed |  4.9s
Computing slice 9/10
[########################################] | 100% Completed |  6.2s
Computing slice 10/10
[########################################] | 100% Completed |  5.6s
{% endraw %} {% raw %}
import moseley as mos 
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax) 
plot_peakslices(x_keVs, slices, ax=ax);

ax1.set_ylim([0, 1.7])

mos.XFluo('Cl', tube_keV=30).plot(ax=ax1, color='b', peak_labels='full')
mos.XFluo('Rh', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
mos.XFluo('K', tube_keV=30).plot(ax=ax1, color='g', peak_labels='simple')
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
titles[n]= '*' + titles[n]
{% endraw %} {% raw %}
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

I would like to inspect the hotmax pixel... We need to check if Rh_L3M5 at 2.701 keV is not the major contributor to this

{% raw %}
ClKa_map = peak_maps[3]
{% endraw %} {% raw %}
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])

y, x, z = hot_pixel

ax.imshow(ClKa_map)
ax1.imshow(imvis_reg)
ax.scatter(x, y, c='r')
ax1.scatter(x, y, c='r');
{% endraw %}

Seems to be a single grain of salt?

How about Lead?

{% raw %}
n = 15

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/8
[########################################] | 100% Completed |  5.8s
Computing slice 1/8
[########################################] | 100% Completed |  5.7s
Computing slice 2/8
[########################################] | 100% Completed |  5.3s
Computing slice 3/8
[########################################] | 100% Completed |  6.5s
Computing slice 4/8
[########################################] | 100% Completed |  5.2s
Computing slice 5/8
[########################################] | 100% Completed | 11.9s
Computing slice 6/8
[########################################] | 100% Completed |  6.0s
Computing slice 7/8
[########################################] | 100% Completed |  5.7s
Computing slice 8/8
[########################################] | 100% Completed |  4.1s
{% endraw %} {% raw %}
import moseley as mos 
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax) 
plot_peakslices(x_keVs, slices, ax=ax);

ax1.set_ylim([0, 1.7])

mos.XFluo('Pb', tube_keV=30).plot(ax=ax1, color='k', peak_labels='full')
mos.XFluo('S', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
{% endraw %} {% raw %}
titles = [f'[{i}]' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

This is most interesting. Let's take a look at the keV map for the overlapping band of sulfur and lead. For sulfur this is the S_KL3 band at 2.311 keV. For lead this is the Pb_M5O2 band and 2.372 keV.

{% raw %}
mid_keV = (2.311 + 2.371) / 2
{% endraw %} {% raw %}
SKa_keV_map = keV_maps[3]
{% endraw %} {% raw %}
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5], constrained_layout=True)

ax.imshow(peak_maps[3])
ax.set_title('Overlapping band peak_maps[3]')
pos = ax1.imshow(SKa_keV_map, vmin=2.2, vmax=2.45)
fig.colorbar(pos, ax=ax1, shrink=0.8)
ax1.set_title('Energy of maximum \nS=2.311 keV vs Pb=2.372 keV');
{% endraw %}

API

Ok, this is the new Gaussian code base. For now we do not return the fancy peak shape parameters to deconvolve overlapping bands. Soon this will be necessary. For now let's simply take a look at the peak maps...

{% raw %}

gaussian[source]

gaussian(x, x0, sigma)

Normal distribution around `x0` with standard deviation `sigma`.
{% endraw %} {% raw %}

fit_gaussian[source]

fit_gaussian(x, y, peak_idx, rel_height=0.2, baseline=None)

Fit single gaussian distribution at `rel_height`.

Returns: `y_gauss`, `baseline`
{% endraw %} {% raw %}

plot_peakslices[source]

plot_peakslices(x, slices, ax=None, y=None, labels=False)

Plot peak slice regions
{% endraw %} {% raw %}

get_peakmaps[source]

get_peakmaps(slices, datastack_file, norm=False)

Integrate peak `slices`  into peak maps and keV maps.

Returns: peak_maps, keV_maps
{% endraw %} {% raw %}

multi_plot[source]

multi_plot(*images, hot_pixel=None, titles=None, roi_list=None, axis_off=False, sharex=True, sharey=True, vmin=None, vmax=None, cmap='viridis', fontsize='medium', zoom_xyc=None, zoom_half_wh=[100, 100])

Inspect multiple images simultaneously...

Fold along multiple rows if n > 4
{% endraw %} {% raw %}

get_slices[source]

get_slices(x_keVs, y_hot, peak_idxs, rel_height=0.2, clip_level=0.05, baseline=None)

Based on single Gaussian peak fit for each peak slice.

Returns `peak_slices`
{% endraw %} {% raw %}

get_hotslices[source]

get_hotslices(datastack_file, n, rel_height=0.2, clip_level=0.05)

Get slices for specifically for hotmax spectrum `n` in `datastack_file`.

Using noiseline as peakbase.


Returns: `peak_slices`
{% endraw %} {% raw %}
{% endraw %}